package jscinttools.scint;
import static java.lang.Math.sqrt;
import static java.lang.Math.pow;
import static java.lang.Math.exp;
import java.lang.System;
import hep.aida.ref.histogram.Measurement;
import hep.aida.IMeasurement;
import hep.aida.IFunctionFactory;
import hep.aida.IFunction;
import hep.aida.ITree;
import jhplot.F1D;
import jscinttools.utilities.aidautilities.AidaFactories;

public class ScintPulseFactory{
	private ScintPulseFactory(){}
	
	public static ScintPulse createScintPulse(double amp, double ampErr, double tau_us, double tauErr_us){
		return createScintPulse(new double[]{amp}, new double[]{ampErr}, new double[]{tau_us}, new double[]{tauErr_us});
	}
	
	public static ScintPulse createScintPulse(double[] amp, double[] ampErr, double[] tau_us, double[] tauErr_us){
		return new ScintPulseImp(amp, ampErr, tau_us, tauErr_us);
	}
}

class ScintPulseImp implements ScintPulse{
	private double[] lyFrac;
	private double[] lyFracErr;
	private double[] amp;
	private double[] ampErr;
	private double[] tau_us;
	private double[] tauErr_us;
	private double lyTot;
	private int size;
	
	public ScintPulseImp(double[] amp, double[] ampErr, double[] tau_us, double[] tauErr_us){
		this.amp = amp;
		this.ampErr = ampErr;
		this.tau_us = tau_us;
		this.tauErr_us = tauErr_us;
		if (((amp.length + ampErr.length + tau_us.length + tauErr_us.length) * 1.0 / 4.0) == amp.length)
			this.size = amp.length;
		calcLYTot();
		lyFracInit();
		//printParam();
	}
	
	private void calcLYTot(){
		for (int i = 0; i < size; i++)
			this.lyTot += amp[i] * tau_us[i];
	}
	
	private void lyFracInit(){
		this.lyFrac = new double[size];
		this.lyFracErr = new double[size];
		for (int i = 0; i < size; i++){
			lyFracComp(i);
		}
	}
	
	private void lyFracComp(int comp){
		this.lyFrac[comp] = amp[comp] * tau_us[comp] / lyTot;
		double sumAiDeltaTaui = calcSumAiDeltaTaui(comp);
		double sumTauiDeltaAi = calcumTauiDeltaAi(comp);
		double t1 = pow(ampErr[comp] * ((1.0 / amp[comp]) + (tau_us[comp] / lyTot)), 2);
		double t2 = pow(tauErr_us[comp] * ((1.0 / tau_us[comp]) + (amp[comp] / lyTot)), 2);
		double t3 = pow(sumAiDeltaTaui / lyTot, 2);
		double t4 = pow(sumTauiDeltaAi / lyTot, 2);
		this.lyFracErr[comp] = this.lyFrac[comp] * sqrt(t1 + t2 + t3 + t4);
	}
	
	private double calcSumAiDeltaTaui(int comp){
		double sumAiDeltaTaui = 0.0;
		for (int i = 0; i < size; i++){
			if (i != comp)
				sumAiDeltaTaui += amp[i] * tauErr_us[i];
		}
		return sumAiDeltaTaui;
	}
	
	private double calcumTauiDeltaAi(int comp){
		double sumTauiDeltaAi = 0.0;
		for (int i = 0; i < size; i++){
			if (i != comp)
				sumTauiDeltaAi += tau_us[i] * ampErr[i];
		}
		return sumTauiDeltaAi;
	}
	
	private void printParam(){
		for (int i = 0; i < size; i++){
			System.out.println("amp_" + i + " = " + amp[i]);
			System.out.println("ampErr_" + i + " = " + ampErr[i]);
			System.out.println("tau_us_" + i + " = " + tau_us[i]);
			System.out.println("tauErr_us_" + i + " = " + tauErr_us[i]);
			System.out.println("lyFrac_" + i + " = " + lyFrac[i]);
			System.out.println("lyFracErr_" + i + " = " + lyFracErr[i]);
		}
	}
	
	public IMeasurement lyFrac(int comp){
		return new Measurement(lyFrac[comp], lyFracErr[comp]);
	}
	
	public IMeasurement lyCorr(double win_us){
		double lyTotTrunc = calcLyTotTrunc(win_us);
		double t1 = sumT1(win_us, lyTotTrunc);
		double t2 = sumT2(win_us, lyTotTrunc);
		double lyCorr = lyTot / lyTotTrunc;
		double lyCorrErr = lyCorr * sqrt(t1 + t2);
		return new Measurement(lyCorr, lyCorrErr);
	}
	
	private double calcLyTotTrunc(double win_us){
		double lyTotTrunc = 0.0;
		for (int i = 0; i < size; i++)
			lyTotTrunc += (amp[i] * tau_us[i] * (1 - exp(-1.0 * win_us / tau_us[i])));
		return lyTotTrunc;
	}
	
	private double sumT1(double win_us, double lyTotTrunc){
		double s = 0.0;
		double t1 = 1.0 / lyTot;
		for (int i = 0; i < size; i++){
			double t2 = (1 - exp(-1.0 * win_us / tau_us[i])) / lyTotTrunc;
			s += pow((t1 - t2) * tau_us[i] * ampErr[i] , 2);
		}
		return s;
	}
	
	private double sumT2(double win_us, double lyTotTrunc){
		double s = 0.0;
		double t1 = 1.0 / lyTot;
		for (int i = 0; i < size; i++){
			double t2 = (1 - (1 + (win_us / tau_us[i])) * exp(-1.0 * win_us / tau_us[i])) / lyTotTrunc;
			s += pow((t1 - t2) * amp[i] * tauErr_us[i], 2);
		}
		return s;
	}
	
	public Coincidence coincEfficiency(double coincWin_us){
		return new CoincidenceImp(coincWin_us);
	}
	
	private class CoincidenceImp implements Coincidence{
		private double coincWin_us;
		private double sumij;
		
		public CoincidenceImp(double coincWin_us){
			this.coincWin_us = coincWin_us;
			this.sumij = sumij();
		}
		
		private double sumij(){
			double s = 0;
			for (int i = 0; i < size; i++){
				for (int j = 0; j < size; j++)
					s += lyFrac[i] * lyFrac[j] * tau_us[j] / (tau_us[i] + tau_us[j]) * exp(-1.0 * coincWin_us / tau_us[j]);
			}
			return 2.0 * s;
		}
		
		public double efficiency(double N1, double N2){
			return 1 - pow(sumij, N1 * N2);
		}
		
		public double efficiency(int N1, int N2){
			return efficiency((double) N1, (double) N2);
		}
		
		public IFunction efficiencyFunc(String pathAndTitle, ITree tree, double chanRatio){
			//String funcDef = "1 - pow(s, x[0])";
			String funcDef = "1 - pow(s, r * pow(x[0] / (r + 1), 2))";
			String param = "s, r";
			String funcName = "Coincidence Efficiency";
			IFunctionFactory ff = AidaFactories.ANALYSIS_FACTORY.createFunctionFactory(tree);
			IFunction func = ff.createFunctionFromScript(pathAndTitle, 1, funcDef, param, funcName);
			func.setParameter("s", sumij);
			func.setParameter("r", chanRatio);
			return func;
		}
		
		public IFunction efficiencyFunc(String title, double chanRatio){
			return efficiencyFunc(title, AidaFactories.TREE, chanRatio);
		}
		
		public F1D efficiencyFunc(String title, double min, double max, double chanRatio){
			return new F1D(efficiencyFunc(title, chanRatio), min, max);
		}
	}	
}
